浅谈 FFT

概述

离散傅里叶变换（Discrete Fourier Transform，缩写为 DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其 DTFT 的频域采样。

FFT 是一种高效实现 DFT 的算法，称为快速傅立叶变换（Fast Fourier Transform，FFT）。它对傅里叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。快速数论变换 (NTT) 是快速傅里叶变换（FFT）在数论基础上的实现。

在 1965 年，Cooley 和 Tukey 发表了快速傅里叶变换算法。事实上 FFT 早在这之前就被发现过了，但是在当时现代计算机并未问世，人们没有意识到 FFT 的重要性。一些调查者认为 FFT 是由 Runge 和 König 在 1924 年发现的。但事实上高斯早在 1805 年就发明了这个算法，但一直没有发表。——OI WIKI

注意：模块并不属于中国计算机学会划定的提高组知识点考察范围。

前置知识

多项式的度

对于一个多项式 $f(x)$，称其最高次项的次数为该多项式的 度（Degree），记作 $\operatorname{deg} f$。

任意角与弧度制

我们在初中学习过角度值，但是角度不是一个数，这给我们深入研究带来了一定的困难，还有其他的问题无法解释清，所以我们换用弧度制描述角。

首先我们用旋转的思路定义角，角可以看成平面内一条射线绕其端点从一个位置旋转到另一个位置形成的图形。开始的位置称为始边，结束的位置称为终边。

我们规定，按 逆时针 方向旋转形成的角叫做 正角，按 顺时针 方向旋转所形成的角叫做 负角，如果这条射线没有做任何旋转，称为 零角。这样我们就把角的概念推向了 任意角。

然后我们介绍 弧度制，把长度等于半径长的弧所对的圆心角称为 1 弧度的角，用符号 $\text{rad}$ 表示，读作：弧度。

一般地，正角的弧度数为正，负角的弧度数为负，零角的弧度数为 0，如果半径为 $r$ 的圆的圆心角 $\alpha$ 所对弧长为 $l$，则 $\vert \alpha \vert =\frac{l}{r} \vert \alpha \vert = \frac{l}{r}$。利用这个公式还可以写出弧长和扇形面积公式，在此略过。

那么，我们发现 $360^\circ$ 的角弧度数为 $2 \pi$，这样有了对应关系之后，我们可以进行角度值和弧度制的转化了。

我们考虑一个角，将其终边再旋转一周，甚至多周，始边位置不动，那么终边位置永远是相同的，我们称这些角为 终边位置相同的角。

与角 $\alpha$ 终边位置相同的角的集合很容易得出，为 $\{ \theta \vert \theta = \alpha + 2k\pi, k \in \mathbb{Z} \}$ 。

可以理解为：给这个角的边不停加转一圈，终边位置不变。

任意角三角函数定义

在平面直角坐标系 $xOy$ 中设 $\angle \beta$ 的始边为 $x$ 轴的正半轴，设点 $P(x, y)$ 为 $\angle \beta$ 的终边上不与原点 $O$ 重合的任意一点，设 $r=OP$，令$\angle \beta = \angle \alpha$，则有 $\sin \alpha = \frac{y}{r}$ ，$\cos \alpha = \frac{x}{r}$ ，$\tan \alpha = \frac{y}{x}$ 。特别地，在 $r = 1$ 时，有 $P(\cos \alpha, \sin \alpha)$ 。

复数

复数的基础知识

由于 $-1$ 不能开方，所以引入了一个数：虚数 $\text{i}$ ，满足 $i^2 = -1$

我们定义形如 $a + b \text{i}$ ，其中 $a,b\in \mathbb{R}$ 的数称为 复数。其中 $\text{i}$ 被称为 虚数单位，全体复数的集合叫做 复数集。

复数通常用 $z$ 表示，即 $z = a + b \text{i}$。这种形式被称为 复数的代数形式。其中 $a$ 称为复数 $z$ 的 实部，表示为 $\operatorname{Re}(z)$ 。$b$ 称为复数 $z$ 的 虚部 ，表示为 $\operatorname{Im}(z)$ 。

我们可以用一个数对 $(a, b)$ 来表示一个复数 $z = a + b \text{i}$ 。可以发现 复数集与平面直角坐标系中的点集一一对应。于是，我们找到了复数的一种几何意义。

那么这个平面直角坐标系就不再一般，因为平面直角坐标系中的点具有了特殊意义——表示一个复数，所以我们把这样的平面直角坐标系称为 复平面，$x$ 轴称为 实轴，$y$ 轴称为 虚轴。我们进一步地说：复数集与复平面内所有的点所构成的集合是一一对应的。

我们考虑到学过的平面向量的知识，发现向量的坐标表示也是一个有序实数对 $(a,b)$，显然，复数 $z=a+b\text{i}$ 对应复平面内的点 $Z(a,b)$，那么它还对应平面向量 $\overrightarrow{OZ}=(a,b)$，于是我们又找到了复数的另一种几何意义：复数集与复平面内的向量所构成的集合是一一对应的（实数 $0$ 与零向量对应）。

于是，我们由向量的知识迁移到复数上来，定义 复数的模 就是复数所对应的向量的模。复数 $z=a+b\text{i}$ 的模 $\vert z \vert = \sqrt{a^2+b^2}$

于是为了方便，我们常把复数 $z=a+b\text{i}$ 称为点 $Z$ 或向量 $\overrightarrow{OZ}$，并规定相等的向量表示同一个复数。

并且由向量的知识我们发现，虚数不可以比较大小（但是实数是可以的）。

复数的三角表示

对于一个复数 $z=a+b\text{i}$ ，我们让 $r = \sqrt{a^2 + b^2}$ ，那么复数 $z$ 可以表示为 $r \cdot ( \cos \theta + \text{i} \sin \theta)$ ，其中 $\cos \theta = \frac{a}{r}$ ，$\sin \theta = \frac{b}{r}$ 。

于是，复数有如下性质（以下 $z_1 = r_1 \cdot (\cos \theta_1 + \text{i} \sin \theta_1)$ ，$z_2 = r_2 \cdot (\cos \theta_2 + \text{i} \sin \theta_2)$）：

1. $z_1 \cdot z_2 = r_1 \cdot r_2 \cdot (\cos (\theta_1 + \theta_2) + \text{i} \sin (\theta_1 + \theta_2))$
2. $z_1^x = z_1 = r_1^x \cdot (\cos (x \cdot \theta_1) + \text{i} \sin (x \cdot \theta_1))$
3. 根据欧拉公式，有 $z_1 = e^{i\pi\theta_1}$
4. $\frac{1}{z_1} = z_1^{-1} = e^{-i\pi\theta_1} = r_1 \cdot (\cos (-\theta_1) + \text{i} \sin (-\theta_1)) = r_1 \cdot (\cos \theta_1 - \text{i} \sin \theta_1)$

单位复根

我们用 $\omega_n$ 表示 $\cos (\frac{2 \pi}{n}) + \text{i} \sin (\frac{2 \pi}{n})$ ，由欧拉公式知 $\omega_n = e^{\frac{2 \pi i}{n}}$ 。

单位复根有三个重要的性质。对于任意正整数 $n$ 和整数 $k$ ，有：

1. $\omega_n^n = 1$
2. $\omega_n^k = \omega_{2n}^{2k}$ ，也就是对于 $2 \vert n$ ，都有 $\omega_n^k = \omega_{\frac n2}^{\frac k2}$
3. $\omega_{2n}^{k + n} = -\omega_{2n}^k$ ，也就是对于 $2 \vert n$ ，都有 $\omega_n^{k + \frac n2} = - \omega_n^k$

这三个性质对后续有重要作用。

多项式的表示

系数表示法

系数表示法就是用一个多项式的各个项系数来表达这个多项式，即使用一个系数序列来表示多项式：

$f(x) = a_0+a_1x+a_2x^2+\cdots +a_{n}x^{n}$

点值表示法

点值表示法是把这个多项式看成一个函数，从上面选取 $n + 1$ 个点，从而利用这 $n + 1$ 个点来唯一地表示这个函数（你可以用高斯消元来理解唯一性）。

$f(x) = a_0+a_1x+a_2x^2+\cdots +a_{n}x^{n} \Leftrightarrow f(x) = \{a_0, a_1, \cdots,a_{n}\}$

设

$\begin{array}{c} f(x_0) = y_0 = a_0 + a_1x_0+a_2x_0^2+a_3x_0^3+ \cdots + a_nx_0^n\\ f(x_1) = y_1 = a_0 + a_1x_1+a_2x_1^2+a_3x_1^3+ \cdots + a_nx_1^n\\ f(x_2) = y_2 = a_0 + a_1x_2+a_2x_2^2+a_3x_2^3+ \cdots + a_nx_2^n\\ \vdots\\ f(x_{n}) = y_{n} = a_0 + a_1x_{n}+a_2x_{n}^2+a_3x_{n}^3+ \cdots + a_nx_{n}^n \end{array}$

那么用点值表示法表示 $f(x)$ 如下：

$f(x) = a_0+a_1x+a_2x^2+\cdots +a_{n}x^{n} \Leftrightarrow f(x) = \{(x_0,y_0),(x_1,y_1), \cdots,(x_n,y_{n})\}$

数论

阶

由欧拉定理可知，对 $a \in \mathbb{Z}$，$m \in \mathbb{N^*}$，若 $\gcd (a,m)=1$，则 $a^{\varphi(m)}\equiv 1 (\operatorname{mod} m)$。

因此满足同余式 $a^n \equiv 1 (\operatorname{mod} m)$ 的最小正整数 $n$ 存在，这个 $n$ 称作 $a$ 模 $m$ 的阶，记作 $\delta_m(a)$。

这里有几个性质：

性质 1：$a,a^2,\cdots,a^{\delta_m(a)}$ 模 $m$ 两两不同余。

性质 2：若 $a^n \equiv 1(\operatorname{mod} m)$，则 $\delta_m(a)∣n$。

性质 3：若 $a^p \equiv a^q (\operatorname{mod} m)$，则有 $p≡q(\operatorname{mod} \delta_m(a))$。

以下两个是与四则运算有关的重要性质。

性质 4：设 $m \in \mathbb{N^*}$，$a,b \in \mathbb{Z}$，$\gcd(a,m)=\gcd(b,m)=1$，则 $\delta_m(ab)=\delta_m(a)\delta_m(b)$ 的充要条件是 $\gcd\big(\delta_m(a),\delta_m(b)\big)=1$

性质 5：设 $k \in \mathbb{N}$，$m \in \mathbb{N^*}$，$a \in \mathbb{Z}$，$\gcd(a,m)=1$，则 $\delta_m(a^k)=\dfrac{\delta_m(a)}{\gcd\big(\delta_m(a),k\big)}$

原根

原根：设 $m \in \mathbb{N^∗}$，$a \in \mathbb{Z}$。若 $\gcd(a,m)=1$，且 $\delta_m(a)=\varphi(m)$，则称 $a$ 为模 $m$ 的原根。特别的，若 $m$ 为质数，则它的原根 $a$ 满足 $a^i \bmod m$ 在 $i \in [1, m - 1]$ 时两两不同。

原根判定定理：若一个数 $g$ 是模 $m$ 的原根，则有对于 $\varphi(m)$ 任何大于 $1$ 且不为自身的因数 $p$，都有 $g^{φ(m)/p} \not\equiv (\operatorname{mod} m)$。

原根个数：若一个数 $m$ 有原根，则它原根的个数为 $\varphi(\varphi(m))$。

原根存在定理：一个数 $m$ 存在原根当且仅当 $m=2,4,p^{\alpha},2p^{\alpha}$，其中 $p$ 为奇素数，$\alpha\in \mathbb{N^*}$。

最小原根的数量级：王元于 1959 年证明了若 $m$ 有原根，其最小原根是不多于 $m^{0.25}$ 级别的。这保证了我们暴力找一个数的最小原根，复杂度是可以接受的。

卷积

对于一个序列，将其中元素一一映射到一个多项式函数的系数上， 这个多项式函数便叫做该序列的生成函数。

形式化地讲，对于序列 $f_{0, 1, \cdots n - 1}$ ，它的生成函数就是 $\sum_{i = 0}^{n - 1} f_i \cdot x^i$ 。

卷积即为生成函数的乘积在对应序列的变换上的的抽象，“卷”即为其作用效果，“积”即为其本质。

对于序列 $f, g$ ，其卷积序列 $f \otimes g$ 满足 $(f \otimes g)_k = \sum_{i = 0}^k f_i \cdot g_{k - i} = \sum_{i + j = k} f_i \cdot g_j$ ，其中 $(f \otimes g)_k$ 表示卷积序列 $f \otimes g$ 的第 $k$ 项。

对于多项式 $f, g$，其多项式的卷积为：$f \otimes g = \sum_{k = 0}^{\operatorname{deg} f + \operatorname{deg} g} ((\sum_{i + j = k} a_i \cdot b_j)\cdot x_k)$ 。

DFT

理论知识

DFT（傅里叶变换）是利用 $O(n \log n)$ 的时间把多项式从系数表示转到了点值表示。

请注意，DFT 只能解决 $n = 2^m$ 时的问题，否则在分治的时候左右不一样长。所以要在第一次 DFT 之前就把序列向上补成长度为 $2^m$（高次系数补 $0$）、最高项次数为 $2^m - 1$ 的多项式。在讲述时，默认 $n = 2^m$ ，其中 $m \in \mathbb{N}$ 。

DFT 它分治地来求当 $x = \omega_n^k$ 的时候 $f(x)$ 的值。也就是

$f(x) = a_0+a_1x+a_2x^2+\cdots +a_{n}x^{n} \Rightarrow f(x) = \{(\omega_n^0,f(\omega_n^0)),(\omega_n^1,f(\omega_n^1)), \cdots,(\omega_n^n,f(\omega_n^n))\}$

它的分治思想体现在将多项式分为奇次项和偶次项处理。举个例子，当 $n = 8$ 时，对于一共 $8$ 项的多项式：

$f(x) = a_0 + a_1x + a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7$

按照次数的奇偶来分成两组，然后右边提出来一个 $x$：

$\begin{aligned} f(x) &= (a_0+a_2x^2+a_4x^4+a_6x^6) + (a_1x+a_3x^3+a_5x^5+a_7x^7)\\ &= (a_0+a_2x^2+a_4x^4+a_6x^6) + x(a_1+a_3x^2+a_5x^4+a_7x^6) \end{aligned}$

分别用奇偶次次项数建立新的函数：

$\begin{aligned} g(x) &= a_0+a_2x+a_4x^2+a_6x^3 \\ h(x) &= a_1+a_3x+a_5x^2+a_7x^3 \end{aligned}$

那么原来的 $f(x)$ 用新函数表示为：

$f(x)=g\left(x^2\right) + x \cdot h\left(x^2\right)$

将 $x = \omega_n^k$ 代入：

$\begin{aligned} f(\omega_n^k) &= g(\omega_n^{2k}) + \omega_n^k \cdot h(\omega_n^{2k}) \\ &= g(\omega_{\frac{n}{2}}^k) + \omega_n^k \cdot h(\omega_{\frac{n}{2}}^k) \end{aligned}$

DFT 最精髓的地方就在于下面这个式子（将 $x = \omega_n^{k + \frac{n}{2}}$ 带入）：

$\begin{aligned} f(\omega_n^{k + \frac{n}{2}}) &= g(\omega_n^{2k + n}) + \omega_n^{k + \frac {n}{2}} \cdot h(\omega_n^{2k + n}) \\ &= g(\omega_n^{2k}) - \omega_n^{k} \cdot h(\omega_n^{2k}) \\ &= g(\omega_{\frac{n}{2}}^k) - \omega_n^{k} \cdot h(\omega_{\frac{n}{2}}^k) \end{aligned}$

也就是说我们只要知道 $g(\omega_{\frac{n}{2}}^k)$ 和 $h(\omega_{\frac{n}{2}}^k)$ 就可以算出 $f(\omega_n^k)$ 和 $f(\omega_n^{k + \frac{n}{2}})$ 了！此时我们只要递归求出 $g(\omega_{\frac{n}{2}}^k)$ 和 $h(\omega_{\frac{n}{2}}^k)$ 即可。

假设算出所有 $f(\omega_n^k)$ 的时间复杂度为 $T(n)$ ，则有 $T(n) = 2 \cdot T(n / 2) + n = O(n \log n)$

在代入值的时候，因为要代入 $n$ 个不同值，所以我们代入 $\omega_n^0, \omega_n^1, \cdots \omega_n^{n - 1}$ 这 $n$ 个不同的值即可。

代码实现方面，STL 提供了复数的模板，当然也可以手动实现。两者区别在于，使用 STL 的 complex 可以调用 exp 函数求出 $\omega_n$。但事实上使用欧拉公式得到的虚数来求 $\omega_n$ 也是等价的。

代码实现

OI WIKI 上的 DFT（STL 版）：

#include <cmath>
#include <complex>

typedef std::complex<double> Comp;  // STL complex

const Comp I(0, 1);  // i
const int MAX_N = 1 << 20;

Comp tmp[MAX_N];

void DFT(Comp *f, int n) {
  if (n == 1) return;
  for (int i = 0; i < n; ++i) tmp[i] = f[i];
  for (int i = 0; i < n; ++i) {  // 偶数放左边，奇数放右边
    if (i & 1)
      f[n / 2 + i / 2] = tmp[i];
    else
      f[i / 2] = tmp[i];
  }
  Comp *g = f, *h = f + n / 2;
  DFT(g, n / 2), DFT(h, n / 2);  // 递归 DFT
  Comp cur(1, 0), step(cos(2 * M_PI / n), sin(2 * M_PI / n));
  // Comp step = exp(I * (2 * M_PI / n)); // 两个 step 定义是等价的
  for (int k = 0; k < n / 2; ++k) { 
    tmp[k] = g[k] + cur * h[k];
    tmp[k + n / 2] = g[k] - cur * h[k];
    cur *= step;
  }
  for (int i = 0; i < n; ++i) f[i] = tmp[i];
}

这里 $g[k] = g(\omega_{\frac{n}{2}}^k)$ ，$h[k] = h(\omega_{\frac{n}{2}}^k)$ ，$step = \omega_n$ ，$cur = \omega_n^k$。

由于 $g(\omega_1^0)$ 和 $h(\omega_1^0)$ 的返回值均为 $a_x$（也就是一个系数），因此我们一开始就可以令 $f[k] = a_k$ ，然后进行 DFT。

IDFT

IDFT（傅里叶逆变换）是利用 $O(n \log n)$ 的时间把多项式从点值表示转到了系数表示。也就是：

$f(x) = \{(\omega_n^0,f(\omega_n^0)),(\omega_n^1,f(\omega_n^1)), \cdots,(\omega_n^n,f(\omega_n^n))\} \Rightarrow f(x) = a_0+a_1x+a_2x^2+\cdots +a_{n}x^{n}$

请注意，这里的 $(\omega_n^k, f(\omega_n^k))$ 与 DFT 中的 $(\omega_n^k, f(\omega_n^k))$ 不同。这里的 $(\omega_n^k, f(\omega_n^k))$ 可以不和 DFT 中的 $(\omega_n^k, f(\omega_n^k))$ 相同。

对此我们有两种理解方式。

理解方式一：线性代数角度

IDFT（傅里叶反变换）的作用，是把目标多项式的点值形式转换成系数形式。而 DFT 本身是个线性变换，可以理解为将目标多项式当作向量，左乘一个矩阵得到变换后的向量，以模拟把单位复根代入多项式的过程：

$\begin{bmatrix}y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix} = \begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}$

我们现在已经知道等号左边 $\begin{bmatrix}y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{bmatrix}$ 的结果了，我们的目标是求出 $\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix}$ 的结果。

根据矩阵的基础知识，我们只要求出中间的大矩阵的逆矩阵即可。

引理：对于一个 $n \times n$ 的矩阵：

$\begin{bmatrix}a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & \cdots & a_{2,n} \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & \cdots & a_{3, n} \\ a_{4, 1} & a_{4, 2} & a_{4, 3} & a_{4, 4} & \cdots & a_{4, n} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & a_{n, 3} & a_{n, 4} & \cdots & a_{n, n} \end{bmatrix}$

满足对于 $\forall i, j \in [1, n]$ ，都有 $a_{i, j} = a_{j, i}$ ，则有

$\begin{bmatrix}a_{1,1} & a_{1,2} & a_{1,3} & a_{1,4} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & \cdots & a_{2,n} \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & \cdots & a_{3, n} \\ a_{4, 1} & a_{4, 2} & a_{4, 3} & a_{4, 4} & \cdots & a_{4, n} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{n, 1} & a_{n, 2} & a_{n, 3} & a_{n, 4} & \cdots & a_{n, n} \end{bmatrix}^{-1} = \frac {1}{n} \begin{bmatrix}a_{1,1}^{-1} & a_{1,2}^{-1} & a_{1,3}^{-1} & a_{1,4}^{-1} & \cdots & a_{1,n}^{-1} \\ a_{2,1}^{-1} & a_{2,2}^{-1} & a_{2,3}^{-1} & a_{2,4}^{-1} & \cdots & a_{2,n}^{-1} \\ a_{3,1}^{-1} & a_{3,2}^{-1} & a_{3,3}^{-1} & a_{3,4}^{-1} & \cdots & a_{3, n}^{-1} \\ a_{4, 1}^{-1} & a_{4, 2}^{-1} & a_{4, 3}^{-1} & a_{4, 4}^{-1} & \cdots & a_{4, n}^{-1} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ a_{n, 1}^{-1} & a_{n, 2}^{-1} & a_{n, 3}^{-1} & a_{n, 4}^{-1} & \cdots & a_{n, n}^{-1} \end{bmatrix}$

而中间这个矩阵 $\begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix}$ 可以用这个引理求得它的逆矩阵。

又 $\omega_n^{-k} = (\omega_n^{-1})^k = \cos (k \cdot \frac{2\pi}{n}) - \text{i} \cdot \sin (k \cdot \frac{2\pi}{n})$

而其它的操作过程与 DFT 是完全相同的。我们可以定义一个函数，在里面加一个参数 $1$ 或者是 $-1$ ，然后把它乘到 $\sin (k \cdot \frac{2\pi}{n})$ 上。传入 $1$ 就是 DFT，传入 $-1$ 就是 IDFT。

理解方式二：单位复根周期性

利用单位复根的周期性同样可以理解 IDFT 与 DFT 之间的关系。

考虑原本的多项式是 $f(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_{n-1} x^{n-1} = \sum_{i = 0}^{n - 1}a_i x^i$。而 IDFT 就是把点值表示还原为系数表示。

考虑 构造法。我们已知 $y_i=f\left( \omega_n^i \right),i\in\{0,1,\cdots,n-1\}$，求 $\{a_0,a_1,\cdots,a_{n-1}\}$。构造多项式如下：

$A(x)=\sum_{i=0}^{n-1}y_ix^i$

相当于把 $\{y_0,y_1,y_2,\cdots,y_{n-1}\}$ 当做多项式 $A$ 的系数表示法。

这时我们有两种推导方式，这对应了两种实现方法。

推导方法一

设 $b_i=\omega_n^{-i}$ ，则多项式 $A$ 在 $x=b_0,b_1,\cdots,b_{n-1}$ 处的点值表示法为 $\left\{ A(b_0),A(b_1),\cdots,A(b_{n-1}) \right\}$ 。

对 $A(x)$ 做一下变换，可以将 $A(b_k)$ 表示为：

$\begin{aligned} A(b_k)&=\sum_{i=0}^{n-1}f(\omega_n^i)\omega_n^{-ik}=\sum_{i=0}^{n-1}\omega_n^{-ik}\sum_{j=0}^{n-1}a_j(\omega_n^i)^{j}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{i(j-k)}=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i \end{aligned}$

下面记 $S(n, a)=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i$

可以发现，在 $a \equiv 0 (\operatorname{mod} n)$ 时，$S(n, a) = n$

而在 $a \not \equiv 0 (\operatorname{mod} n)$ 时，有：

$\begin{aligned} S\left(\omega_n^a\right)=\frac{\left(\omega_n^a\right)^n-\left(\omega_n^a\right)^0}{\omega_n^a-1}=0\end{aligned}$

也就是说，$S(n, a)$ 满足

$S\left(n, a\right)= \left\{\begin{aligned} n &, a \equiv 0 (\operatorname{mod} n) \\ 0 &, a \not \equiv 0 (\operatorname{mod} n) \end{aligned}\right.$

将其带回原式：

$A(b_k)=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(\omega_n^{j-k}\right)^i=\sum_{j=0}^{n-1}a_jS\left(n, j-k\right)=a_k\cdot n$

也就是说给定点 $b_i=\omega_n^{-i}$，则 $A$ 的点值表示法为：

$\begin{aligned} &\left\{ (b_0,A(b_0)),(b_1,A(b_1)),\cdots,(b_{n-1},A(b_{n-1})) \right\}\\ =&\left\{ (b_0,a_0\cdot n),(b_1,a_1\cdot n),\cdots,(b_{n-1},a_{n-1}\cdot n) \right\} \end{aligned}$

说白了，就是说 $\begin{aligned} a_k = \frac{A(\omega_n^{-k})}{n} \end{aligned}$ 。

所以说，我们取单位根为 $\omega_n^{-i}$ ，对着 $A(x)=\sum_{i=0}^{n-1}y_ix^i$ 跑一遍 DFT 后将求得的数都除以 $n$ 即可。

推导方法二

我们直接将 $\omega_n^k$ 代入 $A(x)$。

$\begin{aligned} A(\omega^k_n)&=\sum_{i=0}^{n-1}f(\omega_n^i)\omega_n^{ik}=\sum_{i=0}^{n-1}\omega_n^{ik}\sum_{j=0}^{n-1}a_j(\omega_n^i)^{j}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\omega_n^{i(j+k)}=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(\omega_n^{j+k}\right)^i = A(\omega_n^k) = \sum_{j=0}^{n-1}a_jS(n,j+k) \end{aligned}$

这里 $S(n, a)=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i$

于是，当且仅当 $j + k \equiv 0 (\operatorname{mod} n)$ 时有 $S(n, j + k) = n$ ，其余情况 $S(n, j + k) = 0$ ，因此 $A(\omega_n^k) = a_{n-k}\cdot n$，其中 $k \ge 1$。而 $A(\omega_n^0) = a_0 \cdot n$。

这意味着我们将 $\{y_0,y_1,y_2,\cdots,y_{n-1}\}$ 做 DFT 变换后，$a_1 \sim a_n$ 反转再除以 $n$，同样可以还原 $f(x)$ 的系数表示。

代码实现

现在我们可以集 DFT 和 IDFT 于一身。代码实现如下：

#include <cmath>
#include <complex>

typedef std::complex<double> Comp;  // STL complex

const Comp I(0, 1);  // i
const int MAX_N = 1 << 20;

Comp tmp[MAX_N];

void DFT(Comp *f, int n, int rev) {  // rev = 1, DFT; rev = -1, IDFT
  if (n == 1) return;
  for (int i = 0; i < n; ++i) tmp[i] = f[i];
  for (int i = 0; i < n; ++i) {  // 偶数放左边，奇数放右边
    if (i & 1)
      f[n / 2 + i / 2] = tmp[i];
    else
      f[i / 2] = tmp[i];
  }
  Comp *g = f, *h = f + n / 2;
  DFT(g, n / 2, rev), DFT(h, n / 2, rev);  // 递归 DFT
  Comp cur(1, 0), step(cos(2 * M_PI / n), rev * sin(2 * M_PI / n));
  for (int k = 0; k < n / 2; ++k) {
    tmp[k] = g[k] + cur * h[k];
    tmp[k + n / 2] = g[k] - cur * h[k];
    cur *= step;
  }
  for (int i = 0; i < n; ++i) f[i] = tmp[i];
}

这里 $g[k] = g(\omega_{\frac{n}{2}}^k)$ ，$h[k] = h(\omega_{\frac{n}{2}}^k)$ ，$step = \omega_n$ ，$cur = \omega_n^k$。

由于 $g(\omega_1^0)$ 和 $h(\omega_1^0)$ 的返回值均为 $a_x$（也就是一个系数），因此我们一开始就可以令 $f[k] = a_k$ ，然后进行 DFT。

FFT

位逆序置换

这个算法还可以从“分治”的角度继续优化。我们每一次都会把整个多项式的奇数次项和偶数次项系数分开，一直分到只剩下一个系数。但是，这个递归的过程需要更多的内存。因此，我们可以先“模仿递归”把这些系数在原数组中“拆分”，然后再“倍增”地去合并这些算出来的值。

以 $8$ 项多项式为例，模拟拆分的过程：

• 初始序列为 $\{x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\}$
• 一次拆分之后 $\{x_0, x_2, x_4, x_6\},\{x_1, x_3, x_5, x_7 \}$
• 两次拆分之后 $\{x_0,x_4\} \{x_2, x_6\},\{x_1, x_5\},\{x_3, x_7 \}$
• 三次拆分之后 $\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7 \}$

规律：其实就是原来的那个序列，每个下标用二进制表示，然后把二进制翻转对称一下，就是最终那个位置的下标。比如 $x_1$ 是 001，翻转是 100，也就是 4，而且最后那个位置的下标确实是 4。我们称这个变换为位逆序置换（bit-reversal permutation，国内也称蝴蝶变换）

根据它的定义，我们可以在 $O(n \log n)$ 的时间内求出每个数变换后的结果，以下是代码：（OI WIKI 上的）

void change(Complex y[], int len) {
  int i, j, k;
  for (int i = 1, j = len / 2; i < len - 1; i++) {
    if (i < j) swap(y[i], y[j]);
    // 交换互为小标反转的元素，i<j 保证交换一次
    // i 做正常的 + 1，j 做反转类型的 + 1，始终保持 i 和 j 是反转的
    k = len / 2;
    while (j >= k) {
      j = j - k;
      k = k / 2;
    }
    if (j < k) j += k;
  }
}

我们实际上可以利用 $O(n)$ 递推来实现该变换，以下是代码：（OI WIKI 上的）

// 同样需要保证 len 是 2 的幂
// 记 rev[i] 为 i 翻转后的值
void change(Complex y[], int len) {
  for (int i = 0; i < len; ++i) {
    rev[i] = rev[i >> 1] >> 1;
    if (i & 1) {  // 如果最后一位是 1，则翻转成 len/2
      rev[i] |= len >> 1;
    }
  }
  for (int i = 0; i < len; ++i) {
    if (i < rev[i]) {  // 保证每对数只翻转一次
      swap(y[i], y[rev[i]]);
    }
  }
  return;
}

当然，这个位逆序置换不是没有道理的，可以这么理解：

还是以 $8$ 项多项式为例，模拟拆分的过程（下标用二进制表示）：

• 初始序列为 $\{x_{000}, x_{001}, x_{010}, x_{011}, x_{100}, x_{101}, x_{110}, x_{111}\}$
• 一次拆分之后 $\{x_{000}, x_{010}, x_{100}, x_{110}\},\{x_{001}, x_{011}, x_{101}, x_{111} \}$
• 两次拆分之后 $\{x_{000},x_{100}\} \{x_{010}, x_{110}\},\{x_{001}, x_{101}\},\{x_{011}, x_{111}\}$
• 三次拆分之后 $\{x_{000}\},\{x_{100}\} \{x_{010}\},\{x_{110}\},\{x_{001}\},\{x_{101}\},\{x_{011}\},\{x_{111}\}$

可以发现，每次拆分就是一个按位的基数排序。

第一次拆分之后，第 $0$ 位为 $0$ 的都在左边，为 $1$ 的都在右边。

第二次拆分之后，每一组的第 $1$ 位为 $0$ 的都在左边，为 $1$ 的都在右边。

以此类推，于是就是位逆序置换了。

代码实现（IDFT 推导方法一）

摘自 OI WIKI 的代码：

/*
 * 做 FFT
 * len 必须是 2^k 形式
 * on == 1 时是 DFT，on == -1 时是 IDFT
 */
void fft(Complex y[], int len, int on) {
  change(y, len);
  for (int h = 2; h <= len; h <<= 1) {                  // 模拟合并过程
    Complex wn(cos(2 * PI / h), sin(on * 2 * PI / h));  // 计算当前单位复根 (要根据标记取判断是否取倒数)
    for (int j = 0; j < len; j += h) {
      Complex w(1, 0);  // 计算当前单位复根
      for (int k = j; k < j + h / 2; k++) {
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        y[k] = u + t;  // 这就是把两部分分治的结果加起来
        y[k + h / 2] = u - t;
        // 后半个 “step” 中的ω一定和 “前半个” 中的成相反数
        // “红圈”上的点转一整圈“转回来”，转半圈正好转成相反数
        // 一个数相反数的平方与这个数自身的平方相等
        w = w * wn;
      }
    }
  }
  if (on == -1) {
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
    }
  }
}

代码实现（IDFT 推导方法二）

摘自 OI WIKI 的代码：

/*
 * 做 FFT
 * len 必须是 2^k 形式
 * on == 1 时是 DFT，on == -1 时是 IDFT
 */
void fft(Complex y[], int len, int on) {
  change(y, len);
  for (int h = 2; h <= len; h <<= 1) {             // 模拟合并过程
    Complex wn(cos(2 * PI / h), sin(2 * PI / h));  // 计算当前单位复根 (不用取倒数)
    for (int j = 0; j < len; j += h) {
      Complex w(1, 0);  // 计算当前单位复根
      for (int k = j; k < j + h / 2; k++) {
        Complex u = y[k];
        Complex t = w * y[k + h / 2];
        y[k] = u + t;  // 这就是把两部分分治的结果加起来
        y[k + h / 2] = u - t;
        // 后半个 “step” 中的ω一定和 “前半个” 中的成相反数
        // “红圈”上的点转一整圈“转回来”，转半圈正好转成相反数
        // 一个数相反数的平方与这个数自身的平方相等
        w = w * wn;
      }
    }
  }
  if (on == -1) {
    reverse(y + 1, y + len); // 反转操作
    for (int i = 0; i < len; i++) {
      y[i].x /= len;
    }
  }
}

请注意，如果你的所求序列的数过大，你需要用 long double 才能放置被卡精度（例如 $10^{23}$ 级别）。

“三次变两次”优化

我们发现对于多项式 $P$ 和 $Q$，分别将其放在复数的实部和虚部，有 $(P + Q \text{i})^2 = P^2 - Q^2 + 2(P \cdot Q) \text{i}$。

也就是说，我们可以先将一个多项式的实部和虚部分别放上要卷积的两个多项式 $P$ 和 $Q$ ，再将这个多项式用 FFT，将该多项式平方后利用 IFFT 求出 $(P + Q \text{i})^2$ ，结果就是虚部再除以 $2$ 。

可以发现，整个过程只用了两次 FFT，而前面的都要用三次。这就是“三次变两次”优化。代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1 << 22 | 1;
struct Complex {
	double x, y; // x + yi
	Complex() {x = y = 0.0;}
	Complex(double xx, double yy) {x = xx; y = yy;}
};
Complex operator + (Complex A, Complex B) {
	return Complex(A.x + B.x, A.y + B.y);
}
Complex operator - (Complex A, Complex B) {
	return Complex(A.x - B.x, A.y - B.y);
}
Complex operator * (Complex A, Complex B) {
	return Complex(A.x * B.x - A.y * B.y, A.x * B.y + A.y * B.x);
}
LL A[N], B[N];
Complex fA[N], fB[N], fAns[N], Ans[N];
int Rev[N];
void FFT(Complex *F, int n, int Sign) { // Sign = 1: FFT , Sign = -1: IFFT
	for (int i = 0; i < n; i++) {
		if (i < Rev[i]) swap(F[i], F[Rev[i]]);
	}
	for (int len = 2; len <= n; len <<= 1) {
		int mid = len >> 1;
		Complex wl = Complex(cos(M_PI / mid), Sign * sin(M_PI / mid));
		for (int i = 0; i < n; i += len) {
			Complex Ang = Complex(1, 0);
			for (int j = 0; j < mid; j++) {
				Complex tx = F[i + j], ty = F[i + j + mid];
				F[i + j] = tx + Ang * ty; F[i + j + mid] = tx - Ang * ty;
				Ang = Ang * wl;
			}
		}
	}
}
int main() {
	int n, m, p = 0;
	scanf("%d%d", &n, &m);
	for (int i = 1; i <= n + 1; i++) scanf("%lld", A + i), fA[i - 1].x = A[i];
	for (int i = 1; i <= m + 1; i++) scanf("%lld", B + i), fA[i - 1].y = B[i];
	while ((1 << p) <= n + m) p++;
	for (int i = 1, End = (1 << p); i < End; i++) {
		Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (p - 1));
	}
	FFT(fA, 1 << p, 1);
	for (int i = 0, End = (1 << p); i < End; i++) {
		fAns[i] = fA[i] * fA[i];
	}
	FFT(fAns, 1 << p, -1);
	for (int i = 0; i <= n + m; i++) {
		printf("%d%c", (int)(fAns[i].y / (1 << p) / 2.0 + 0.5), " \n"[i == n + m]); // 四舍五入防止卡精度
	}
	return 0;
}

NTT

NTT（快速数论变换）解决的是多项式乘法带模数的情况，可以说有些受模数的限制，数也比较大。但是比较方便，毕竟没有复数部分。

由于 FFT 中有浮点数乘法，常数较大。而 NTT 中均为整数，常数相对较小。

理论知识

对于一个质数 $p$，若它满足 $p = q \cdot n + 1$ ，其中 $n$ 满足 $n = 2^m,\ m \in \mathbb{N}$，那么它的原根 $g^{qn} \equiv 1 (\operatorname{mod} p)$。

我们可以用 $g_i$ 表示 $g^{\frac{p - 1}{i}}$ 。将 $g_i$ 看作 $\omega_i$ 的等价，可以发现 $g_i$ 满足与 $\omega_i$ 相似的性质，例如对于任意正整数 $n = 2^m$ 和整数 $k$ ，有：

1. $\omega_n^n = 1 \Leftrightarrow g_n^n \equiv 1 (\operatorname{mod} p)$

2. $\omega_n^k = \omega_{2n}^{2k} \Leftrightarrow g_n^k \equiv g_{2n}^{2k} (\operatorname{mod} p)$

3. $\omega_{2n}^{k + n} = -\omega_{2n}^k \Leftrightarrow g_{2n}^{k + n} \equiv -g_{2n}^{k} (\operatorname{mod} p)$

4. $g_n^{\frac{n}{2}} \equiv -1 (\operatorname{mod} p)$

这里常见的模数 $p$ 有：

$p=998244353=7 \times 17 \times 2^{23}+1, g=3$

$p = 1004535809 = 479 \times 2^{21}+1, g=3$

下面这个也能用：

$p = 469762049 = 7 \times 2^{26} + 1, g = 3$

这里 $g$ 表示 $p$ 的一个原根。

NTT 部分与 FFT 一样，这里推导一下 INTT。

INTT 推导

和 IDFT 一样，考虑 构造法。我们已知 $y_i=f\left(g_n^i \right),i\in\{0,1,\cdots,n-1\}$，求 $\{a_0,a_1,\cdots,a_{n-1}\}$。构造多项式如下：

$A(x)=\sum_{i=0}^{n-1}y_ix^i$

相当于把 $\{y_0,y_1,y_2,\cdots,y_{n-1}\}$ 当做多项式 $A$ 的系数表示法。

我们直接将 $g_n^k$ 代入 $A(x)$。

$\begin{aligned} A(g^k_n)&=\sum_{i=0}^{n-1}f(g_n^i)g_n^{ik}=\sum_{i=0}^{n-1}g_n^{ik}\sum_{j=0}^{n-1}a_j(g_n^i)^{j}\\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_jg_n^{i(j+k)}=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(g_n^{j+k}\right)^i = A(g_n^k) = \sum_{j=0}^{n-1}a_jS(n,j+k) \end{aligned}$

这里 $S(n, a)=\sum_{i=0}^{n-1}\left(g_n^a\right)^i$

在 $a \equiv 0 (\operatorname{mod} n)$ 时，$S(n, a) = n \equiv n (\operatorname{mod} p)$ 。

在 $a \not \equiv 0 (\operatorname{mod} n)$ 时，$\begin{aligned} S(n, a) = \frac{(g_n^a)^n - (g_n^a)^0}{g_n^a - 1} \equiv \frac{(g_n^n)^a - 1}{g_n^a - 1} \equiv 0 (\operatorname{mod} p) \end{aligned}$ 。

于是，当且仅当 $j + k \equiv 0 (\operatorname{mod} n)$ 时有 $S(n, j + k) \equiv n (\operatorname{mod} p)$ ，其余情况 $S(n, j + k) \equiv 0 (\operatorname{mod} p)$ ，因此 $A(g_n^k) \equiv a_{n-k}\cdot n (\operatorname{mod} p)$，其中 $k \ge 1$。而 $A(g_n^0) \equiv a_0 \cdot n (\operatorname{mod} p)$ 。

这意味着我们将 $\{y_0,y_1,y_2,\cdots,y_{n-1}\}$ 做 NTT 变换后，反转 $a_1 \sim a_n$ 后再乘以 $n$ 的逆元，同样可以还原 $f(x)$ 的系数在模 $p$ 后的表示。

代码实现

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1 << 22 | 1;
const LL mod = 998244353;
LL A[N], B[N];
LL fA[N], fB[N], fAns[N], Ans[N];
int Rev[N];
inline LL qpow(LL n, LL p, LL mod = 998244353ll) { // 求逆元用
    LL ret = 1ll;
    while (p) {
        if (p & 1) (ret *= n) %= mod;
        (n *= n) %= mod; p >>= 1;
    }
    return ret;
}
void NTT(LL *F, int n, int Sign) { // Sign = 1: NTT , Sign = -1: INTT
	for (int i = 0; i < n; i++) {
		if (i < Rev[i]) swap(F[i], F[Rev[i]]);
	}
	for (int len = 2; len <= n; len <<= 1) {
		int mid = len >> 1;
        LL wl = qpow(3, (mod - 1) / len);
		for (int i = 0; i < n; i += len) {
			LL Ang = 1;
			for (int j = 0; j < mid; j++) {
				LL tx = F[i + j], ty = (F[i + j + mid] * Ang) % mod;
				F[i + j] = (tx + ty) % mod; F[i + j + mid] = (tx - ty + mod) % mod;
				Ang = (Ang * wl) % mod;
			}
		}
	}
    if (Sign == -1) {
        for (int i = 1; i < (n >> 1); i++) swap(F[i], F[n - i]);
        LL Invn = qpow(n, mod - 2);
        for (int i = 0; i < n; i++) {
            (F[i] *= Invn) %= mod;
        }
    }
}
int main() {
	int n, m, p = 0;
	scanf("%d%d", &n, &m);
	for (int i = 1; i <= n + 1; i++) scanf("%lld", A + i), fA[i - 1] = A[i];
	for (int i = 1; i <= m + 1; i++) scanf("%lld", B + i), fB[i - 1] = B[i];
	while ((1 << p) <= n + m) p++;
	for (int i = 1, End = (1 << p); i < End; i++) {
		Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (p - 1));
	}
	NTT(fA, 1 << p, 1); NTT(fB, 1 << p, 1);
	for (int i = 0, End = (1 << p); i < End; i++) {
		fAns[i] = (fA[i] * fB[i]) % mod;
	}
	NTT(fAns, 1 << p, -1);
    for (int i = 0; i <= n + m; i++) {
        printf("%lld%c", fAns[i], " \n"[i == n + m]);
    }
	return 0;
}

可以发现，用 NTT 交模板题比用 FFT 交模板题的评测用时短，内存用的少，也少去了手写复数的麻烦，是写多项式的不错的选择。

例题

【模板】多项式乘法（FFT）

【模板】A*B Problem升级版（FFT快速傅里叶）

[ZJOI2014]力

题目链接

给定 $n$ 个数 $q_1, q_2, \cdots q_n$ ，定义 $\begin{aligned} F_j = \sum_{i = 1}^{j - 1} \frac{q_i \cdot q_j}{(j - i) ^ 2} - \sum_{i = j + 1}^n \frac{q_i \cdot q_j}{(j - i)^2} \end{aligned}$ ，以及 $\begin{aligned} E_i = \frac{F_i}{q_i} \end{aligned}$ 。现在请你对 $1 \le i \le n$ ，求 $E_i$ 的值。

$1 \le n \le 10^5$ ，$0 < q_i < 10^9$ 。

样例输入：3 3.14 19.26 8.17 ，输出：-21.302 -5.030 20.045 。

题目解答

这题是一个凑卷积的题。

可以发现，这题中 $\begin{aligned} E_i = \sum_{j = 1}^{i - 1} \frac{q_j}{(j - i) ^ 2} - \sum_{j = i + 1}^n \frac{q_j}{(j - i)^2} \end{aligned}$ 。我们将 $E_i$ 拆成两部分：$\begin{aligned} A_i = \sum_{j = 1}^{i - 1} \frac{q_j}{(j - i) ^ 2} \end{aligned}$ 和 $\begin{aligned} B_i = \sum_{j = i + 1}^n \frac{q_j}{(j - i)^2} \end{aligned}$ 。

对于 $A_i$ ，有$\begin{aligned} A_i = \sum_{j = 1}^{i - 1} \frac{q_j}{(j - i) ^ 2} = \sum_{j = 1}^{i - 1} q_j \cdot \frac{1}{(i - j)^2} \end{aligned}$ 。

定义 $\begin{aligned} g_i = \frac{1}{i^2} \end{aligned}$ ，那么有 $\begin{aligned} A_i = \sum_{j = 1}^{i - 1} q_j \cdot \frac{1}{(i - j)^2} = \sum_{j = 1}^{i - 1} q_j \cdot g_{i - j} \end{aligned}$ ，是个卷积形式，可以用 FFT 进行加速。

对于 $B_i$ ，有 $\begin{aligned} B_i = \sum_{j = i + 1}^n \frac{q_j}{(j - i)^2} = \sum_{j = i + 1}^n q_j \cdot \frac{1}{(j - i)^2} \end{aligned}$

还是定义 $\begin{aligned} g_i = \frac{1}{i^2} \end{aligned}$ ，那么有 $\begin{aligned} B_i = \sum_{j = i + 1}^{n} q_j \cdot \frac{1}{(j - i)^2} = \sum_{j = i + 1}^{n} q_j \cdot g_{j - i} \end{aligned}$

可以发现，这并不是一个标准的卷积形式。我们可以将序列 $q$ 进行反转。令 $p_i = q_{n + 1 - i}$ 。那么有 $\begin{aligned} g_i = \frac{1}{i^2} \end{aligned}$ ，那么有 $\begin{aligned} B_i = \sum_{j = 1}^{i - 1} q_j \cdot g_{j - i} = \sum_{j = i + 1}^{n} p_{n + 1 - j} \cdot g_{j - i} \end{aligned}$

乍一看这不是卷积的形式，不过如果令 $a = n + 1 - j$ ，那么 $\begin{aligned} B_i = \sum_{j = i + 1}^{n} p_{n + 1 - j} \cdot g_{j - i} = \sum_{a = 1}^{n - i} p_a \cdot g_{(n + 1 - i) - a} \end{aligned}$ 。

于是这又是一个卷积形式，用 FFT 就做完了。

最后将 $A_i$ 和 $B_i$ 分别计算后，$E_i = A_i - B_{n + 1 - i}$ 。

[SDOI2015]序列统计

题目链接

给定一个集合 $S$ ，满足 $S \subseteq \{0, 1, \cdots, m - 1\}$ ，求从这个集合 $S$ 中有顺序地选取 $n$ 个数（可以重复），使得这 $n$ 个数的乘积 $\operatorname{mod} m$ 为 $x$ 的方案数 模上 $1004535809$ 的值。其中 $m$ 为质数。

对于 $100\%$ 的数据，$1 \le n \le 10^9$，$3 \le m \le 8000$，$1 \le x < m$。保证 $m$ 为质数，且 $S$ 中的元素两两不同。

样例输入：4 3 1 2\n1 2 ，输出：8

样例解释：可以生成的满足要求的不同的数列有(1,1,1,1)、(1,1,2,2)、(1,2,1,2)、(1,2,2,1)、(2,1,1,2)、(2,1,2,1)、(2,2,1,1)、(2,2,2,2)。

题目解答

这是个较为经典的问题。下面引入一些简单的题目来进行解答。

例题 1

给定一个集合 $S$，满足 $S \subseteq \{ 0, 1, \cdots , n \}$ ，在这个集合里选取 $3$ 个数（可以重复），使得加和为 $A$ 的方案数。

例题 1 解答

可以运用桶的思想。设非负整数 $x$ 在集合 $S$ 中出现了 $h(x)$ 次，那么可以得到答案为 $\begin{aligned} \sum_{i = 0}^A \sum_{j = 0}^{A - i} h(i) \times h(j) \times h(A - i - j) \end{aligned}$ 。

时间复杂度是 $O(n^2)$ ，有点慢。不过我们可以先求出从 $S$ 中选取 $2$ 个数的方案。设从 $S$ 中选取 $2$ 个数的加和为 $x$ 的方案数为 $h_2(x)$ ，则可以得到答案为 $\begin{aligned} \sum_{i = 0}^A h(i) \times h_2(A - i) \end{aligned}$ 。可以 $O(n)$ 求出。

那么如何快速求出 $h_2$ 呢？可以发现，$\begin{aligned} h_2(x) = \sum_{i = 0}^x h(i) \times h_2(x - i) \end{aligned}$ 。这是个标准的卷积形式，可以使用 FFT/NTT 快速求出 $h_2$ 。

这样，总的复杂度就是 $O (n \log n)$ 了。

例题 2

给定一个集合 $S$，满足 $S \subseteq \{ 0, 1, \cdots , n \}$ ，在这个集合里选取 $k$ 个数（可以重复），使得加和为 $A$ 的方案数。

例题 2 解答

类比于例题 1，不难得到一种做法：设从 $S$ 中选取 $t$ 个数的加和为 $x$ 的方案数为 $h_t(x)$ ，则有这样的递推式：$\begin{aligned} h_t(x) = \sum_{i = 0}^A h(i) \times h_{t - 1}(x - i) \end{aligned}$

可以发现，这样做的话复杂度为 $O(k A \log A)$ 的，还可以在进行优化。

可以考虑倍增。可以发现，选取 $2^r$ 个数相当于在选取 $2^{r - 1}$ 个数的 可能的和 的集合中选取 $2$ 个数。那么就有递推式：$\begin{aligned} h_{2^r}(x) = \sum_{i = 0}^A h_{2^{r - 1}}(i) \times h_{2^{r - 1}}(x - i) \end{aligned}$

如此这般，就可以预处理处所有的 $h_{2^r}$ ，之后进行与快速幂类似的操作，将 $k$ 二进制拆解，可以在 $O(A \log A \log k)$ 的复杂度下实现。

CF755G PolandBall and Many Other Balls

题目链接

有 $n$ 个球在一排上，定义一个组可以只包含 $1$ 个球或者包含 $2$ 个相邻的球，而一个球最多可以分到一个组中，对于 $\forall m \in [1, k]$ ，求从这些球中取出 $m$ 组的方案数。

$1 \le n \le 10^5$ ，$1 \le k < 2^{15}$ 。

样例输入：3 3 ，输出：5 5 1

样例解释：
$k = 1$ 时，有 $(1), (2), (3), (12), (23)$ ，共 $5$ 种方案。
$k = 2$ 时，有 $(12)(3),(1)(23),(1)(2),(1)(3),(2)(3)$ ，共 $5$ 种方案。
$k = 3$ 时，有 $(1)(2)(3)$ ，共 $1$ 种方案。

题目解答

这是个利用 FFT/NTT 来优化 dp 的题目。

拓展

多项式求逆

题目链接

给定一个多项式 $f(x)$ ，请求出一个多项式 $g(x)$， 满足 $f(x) \cdot g(x) \equiv 1 (\operatorname{mod} x^n)$。系数对 $998244353$ 取模。其中 $\operatorname{deg} f = \operatorname{deg} g = n - 1$。

这里，若一个多项式 $f(x) \equiv 1 (\operatorname{mod} x^n)$ ，则 $f(x) = g(x) \cdot x^n + 1$ ，其中 $g(x)$ 为关于 $x$ 的多项式（这一行的 $f, g$ 与题干中的无关）。

题目解答

利用倍增法。加速递推就可以用倍增。

可以发现，当 $f(x) = c$ 时，$g(x) \equiv c^{-1}(\operatorname{mod} 998244353)$ 。

假设现在已经求出了 $f(x)$ 在模 $x^{\left\lceil\frac{n}{2}\right\rceil}$ 意义下的逆元 $h(x)$ ，要求 $g(x)$，则有：

$f(x) \cdot h(x) \equiv 1 (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil}) \cdots (8.1)$

又 $f(x) \cdot g(x) \equiv 1 (\operatorname{mod} x^n) \Rightarrow f(x) \cdot g(x) \equiv 1 (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil}) \cdots (8.2)$

将两式 $(8.1)$ 与 $(8.2)$ 相减，得到 $f(x) \cdot (h(x) - g(x)) \equiv 0 (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil}) \cdots (8.3)$

由于 $f(x) \cdot h(x) \equiv 1 (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil})$ 知 $f(x) \not \equiv 0 (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil})$

所以可由 $(8.3)$ 得出 $h(x) - g(x) \equiv 0 (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil})$

将其两边平方，得到 $(h(x) - g(x))^2 \equiv 0 (\operatorname{mod} x^n)$

将左边展开，得到 $(h(x))^2 - 2 \cdot g(x) \cdot h(x) + (g(x))^2 \equiv 0 (\operatorname{mod} x^n)$

移项，得到 $(g(x))^2 \equiv 2 \cdot g(x) \cdot h(x) - (h(x))^2 (\operatorname{mod} x^n)$

两边同时除以一个 $g(x)$，得到 $g(x) \equiv 2 \cdot h(x) - (h(x))^2 \cdot (g(x))^{-1} (\operatorname{mod} x^n)$。

发现 $(g(x))^{-1} \equiv f(x) (\operatorname{mod} x^n)$

所以 $g(x) \equiv 2 \cdot h(x) - (h(x))^2 \cdot f(x) = h(x) \cdot (2 - h(x) \cdot f(x)) (\operatorname{mod} x^n)$

用 NTT 即可求出 $g(x)$ 。

时间复杂度 $T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right)$ 。

这里贴一下代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1 << 22 | 5;
const LL mod = 998244353ll;
LL A[N], B[N]; // f, g
LL fA[N], fB[N], fAns[N];
int rev[23][N];
inline LL qpow(LL n, LL p, LL mod = 998244353ll) {
    LL ret = 1ll;
    while (p) {
        if (p & 1) (ret *= n) %= mod;
        (n *= n) %= mod; p >>= 1;
    }
    return ret;
}
void NTT(LL *F, int *Rev, int n, int Sign) { // Sign = 1: NTT , Sign = -1: INTT
	for (int i = 0; i < n; i++) {
		if (i < Rev[i]) swap(F[i], F[Rev[i]]);
	}
	for (int len = 2; len <= n; len <<= 1) {
		int mid = len >> 1;
        LL wl = qpow(3, (mod - 1) / len);
		for (int i = 0; i < n; i += len) {
			LL Ang = 1;
			for (int j = 0; j < mid; j++) {
				LL tx = F[i + j], ty = (F[i + j + mid] * Ang) % mod;
				F[i + j] = (tx + ty) % mod; F[i + j + mid] = (tx - ty + mod) % mod;
				Ang = (Ang * wl) % mod;
			}
		}
	}
    if (Sign == -1) {
        for (int i = 1; i < (n >> 1); i++) swap(F[i], F[n - i]);
        LL Invn = qpow(n, mod - 2);
        for (int i = 0; i < n; i++) {
            (F[i] *= Invn) %= mod;
        }
    }
}
LL Tmp[N]; // Tmp 用来存 f
void Polyinv(LL *f, LL *g, int n) {
    if (n == 1) {g[0] = qpow(f[0], mod - 2); return;}
    Polyinv(f, g, (n + 1) >> 1);
    int p = 0; while ((1 << p) <= n + n) p++;
    for (int i = 0; i < (1 << p); i++) Tmp[i] = (i < n) ? f[i] : 0; // 超过 n 的都设为 0
    NTT(Tmp, rev[p], 1 << p, 1); NTT(g, rev[p], 1 << p, 1);
    for (int i = 0; i < (1 << p); i++) Tmp[i] = ((g[i] * ((2ll - Tmp[i] * g[i]) % mod)) % mod + mod) % mod;
    NTT(Tmp, rev[p], 1 << p, -1);
    for (int i = 0; i < (1 << p); i++) g[i] = (i < n) ? Tmp[i] : 0; // 超过 n 的都设为 0, 否则可能影响后面结果.
}
int main() {
	int n; scanf("%d", &n);
	for (int i = 0; i < n; i++) scanf("%lld", A + i);
    for (int k = 1; k < 22; k++) {
        for (int i = 1, End = (1 << k); i < End; i++) {
            rev[k][i] = (rev[k][i >> 1] >> 1) | ((i & 1) << (k - 1));
        }
    }
    Polyinv(A, B, n);
    for (int i = 0; i < n; i++) {
        printf("%lld%c", B[i], " \n"[i == n - 1]);
    }
	return 0;
}

多项式开根

题目链接

给定一个多项式 $f(x)$ ，请求出一个多项式 $g(x)$， 满足 $(g(x))^2 \equiv f(x) (\operatorname{mod} x^n)$。系数对 $998244353$ 取模。其中 $\operatorname{deg} f = \operatorname{deg} g = n - 1$。该题保证多项式的 $0$ 次项为 $1$ （否则需要用到二次剩余）。

题目解答

还是利用倍增法。可以发现，当 $f(x) = 1$ 时，$g(x) = 1$ 。

假设现在已经求出了 $f(x)$ 在模 $x^{\left\lceil\frac{n}{2}\right\rceil}$ 意义下的平方根 $h(x)$ ，要求 $g(x)$，则有：

$(h(x))^2 \equiv f(x) (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil}) \Rightarrow (h(x))^2 - f(x) \equiv 0 (\operatorname{mod} x^{\left\lceil\frac{n}{2}\right\rceil})$

将两边平方，得到 $((h(x))^2 - f(x))^2 \equiv 0 (\operatorname{mod} x^n)$

通过一定的变换可以得到 $((h(x))^2 + f(x))^2 \equiv 4 \cdot (h(x))^2 \cdot f(x) (\operatorname{mod} x^n)$ 。

于是 $\left(\frac{(h(x))^2+f(x)}{2 \cdot h(x)}\right)^{2} \equiv f(x) (\operatorname{mod} x^n)$

两边同时开根，得到 $\left(\frac{(h(x))^2+f(x)}{2 \cdot h(x)}\right)^{2} \equiv g(x) (\operatorname{mod} x^n)$

所以说 $g(x) \equiv 2^{-1} \cdot h(x) + 2^{-1} \cdot h^{-1}(x) \cdot f(x) (\operatorname{mod} x^n)$ 。

用 NTT 即可求出 $g(x)$ 。

时间复杂度 $T\left(n\right)=T\left(\frac{n}{2}\right)+O\left(n\log{n}\right)=O\left(n\log{n}\right)$ 。

这里贴一下代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1 << 22 | 5;
const LL mod = 998244353ll, inv2 = 499122177ll;
LL A[N], B[N]; // f, g
LL fA[N], fB[N], fAns[N];
int rev[23][N];
inline LL qpow(LL n, LL p, LL mod = 998244353ll) {
    LL ret = 1ll;
    while (p) {
        if (p & 1) (ret *= n) %= mod;
        (n *= n) %= mod; p >>= 1;
    }
    return ret;
}
void NTT(LL *F, int *Rev, int n, int Sign) { // Sign = 1: NTT , Sign = -1: INTT
	for (int i = 0; i < n; i++) {
		if (i < Rev[i]) swap(F[i], F[Rev[i]]);
	}
	for (int len = 2; len <= n; len <<= 1) {
		int mid = len >> 1;
        LL wl = qpow(3, (mod - 1) / len);
		for (int i = 0; i < n; i += len) {
			LL Ang = 1;
			for (int j = 0; j < mid; j++) {
				LL tx = F[i + j], ty = (F[i + j + mid] * Ang) % mod;
				F[i + j] = (tx + ty) % mod; F[i + j + mid] = (tx - ty + mod) % mod;
				Ang = (Ang * wl) % mod;
			}
		}
	}
    if (Sign == -1) {
        for (int i = 1; i < (n >> 1); i++) swap(F[i], F[n - i]);
        LL Invn = qpow(n, mod - 2);
        for (int i = 0; i < n; i++) {
            (F[i] *= Invn) %= mod;
        }
    }
}
LL Tmp[N]; // 用来存 f
void Polyinv(LL *f, LL *g, int n) {
    if (n == 1) {g[0] = qpow(f[0], mod - 2); return;}
    Polyinv(f, g, (n + 1) >> 1);
    int p = 0; while ((1 << p) <= n + n) p++;
    for (int i = 0; i < (1 << p); i++) Tmp[i] = (i < n) ? f[i] : 0;
    NTT(Tmp, rev[p], 1 << p, 1); NTT(g, rev[p], 1 << p, 1);
    for (int i = 0; i < (1 << p); i++) Tmp[i] = ((g[i] * ((2ll - Tmp[i] * g[i]) % mod)) % mod + mod) % mod;
    NTT(Tmp, rev[p], 1 << p, -1);
    for (int i = 0; i < (1 << p); i++) g[i] = (i < n) ? Tmp[i] : 0;
}
LL ginv[N]; // 用来存 h 的逆元
void PolySqrt(LL *f, LL *g, int n) {
    if (n == 1) {g[0] = 1; return;}
    PolySqrt(f, g, (n + 1) >> 1);
    int p = 0; while ((1 << p) <= n + n) p++;
    for (int i = 0; i < (1 << p); i++) ginv[i] = 0;
    Polyinv(g, ginv, n);
    for (int i = 0; i < (1 << p); i++) Tmp[i] = (i < n) ? f[i] : 0;
    NTT(Tmp, rev[p], 1 << p, 1); NTT(g, rev[p], 1 << p, 1); NTT(ginv, rev[p], 1 << p, 1);
    for (int i = 0; i < (1 << p); i++) Tmp[i] = (inv2 * ((g[i] + ginv[i] * Tmp[i]) % mod)) % mod;
    NTT(Tmp, rev[p], 1 << p, -1);
    for (int i = 0; i < (1 << p); i++) g[i] = (i < n) ? Tmp[i] : 0;
}
int main() {
	int n; scanf("%d", &n);
	for (int i = 0; i < n; i++) scanf("%lld", A + i);
    for (int k = 1; k < 22; k++) {
        for (int i = 1, End = (1 << k); i < End; i++) {
            rev[k][i] = (rev[k][i >> 1] >> 1) | ((i & 1) << (k - 1));
        }
    }
    PolySqrt(A, B, n);
    for (int i = 0; i < n; i++) {
        printf("%lld%c", B[i], " \n"[i == n - 1]);
    }
	return 0;
}

请注意，如果该题多项式的 $0$ 次项不为 $1$ ，就需要使用二次剩余来计算 $g(x)$ 的 $0$ 次项。

任意模数多项式乘法 (MTT)

题目链接

题目大意

给定两个多项式 $f(x), g(x)$ ，请求出 $f(x) * g(x)$ 。系数对 $p$ 取模，且不保证 $p$ 可以分解成 $p = a \cdot 2^k + 1$ 之形式。

三模数 NTT

任意模数 NTT，最大的数为 $p^2\times\max\{n,m\}\leq 10^{23}$ ，所以一般选三个模数即可，求出这三个模数下的答案，然后中国剩余定理即可。

假设这一位需要求得的答案是 $x$，三个模数分别为 $A,B,C$ ，那么：

$x \equiv x_1 (\operatorname{mod} A),\ x \equiv x_2 (\operatorname{mod} B),\ x \equiv x_3 (\operatorname{mod} C)$

先解出前两个：$x \equiv x_1 (\operatorname{mod} A),\ x \equiv x_2 (\operatorname{mod} B) \Rightarrow$ 设 $x_1 + k_1 \cdot A = x_2 + k_2 \cdot B$ 有 $x_1 + k_1 \cdot A \equiv x_2 (\operatorname{mod} B)$。

所以 $k_1 \equiv \frac{x_2 - x_1}{A} (\operatorname{mod} B)$，这样就求出了 $k_1$。也就求出了 $x \equiv x_1 + k_1 \cdot A (\operatorname{mod} AB)$，令 $x_4=x_1+k_1 \cdot A$。

现在就变成了以下两个式子：$x \equiv x_4 = x_1 + k_1 \cdot A (\operatorname{mod} AB),\ x \equiv x_3 (\operatorname{mod} C)$。

设 $x_4 + k_4 \cdot AB = x_3 + k_3 \cdot C$ 有 $x_4 + k_4 \cdot AB \equiv x_3 (\operatorname{mod} C)$。

所以 $k_4 \equiv \frac{x_3 - x_4}{AB} (\operatorname{mod} C)$，这样就求出了 $k_4$。也就求出了 $x \equiv x_4 + k_4 \cdot AB (\operatorname{mod} ABC)$。因为我们要求的 $x$ 小于 $ABC$，而 $x_4 + k_4 \cdot AB < ABC$，所以 $x = x_4 + k_4 \cdot AB$，就做完了。共 9 次 NTT，常数巨大，最慢的点 847ms。代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1 << 22 | 5;
const LL mod1 = 998244353ll, mod2 = 1004535809ll, mod3 = 469762049ll;
LL A[N], B[N], mod;
int rev[23][N];
inline LL qpow(LL n, LL p, LL mod = 998244353ll) {
    LL ret = 1ll;
    while (p) {
        if (p & 1) (ret *= n) %= mod;
        (n *= n) %= mod; p >>= 1;
    }
    return ret;
}
void NTT(LL *F, int *Rev, int n, int Sign, LL mod = 998244353ll) { // Sign = 1: NTT , Sign = -1: INTT
	for (int i = 0; i < n; i++) {
		if (i < Rev[i]) swap(F[i], F[Rev[i]]);
	}
	for (int len = 2; len <= n; len <<= 1) {
		int mid = len >> 1;
        LL wl = qpow(3, (mod - 1) / len, mod);
		for (int i = 0; i < n; i += len) {
			LL Ang = 1;
			for (int j = 0; j < mid; j++) {
				LL tx = F[i + j], ty = (F[i + j + mid] * Ang) % mod;
				F[i + j] = (tx + ty) % mod; F[i + j + mid] = (tx - ty + mod) % mod;
				Ang = (Ang * wl) % mod;
			}
		}
	}
    if (Sign == -1) {
        for (int i = 1; i < (n >> 1); i++) swap(F[i], F[n - i]);
        LL Invn = qpow(n, mod - 2, mod);
        for (int i = 0; i < n; i++) {
            (F[i] *= Invn) %= mod;
        }
    }
}
LL A1[N], A2[N], A3[N];
LL B1[N], B2[N], B3[N];
LL C1[N], C2[N], C3[N];
LL F(LL c1, LL c2, LL c3, LL mod) {
    // mod1, mod2, mod3, mod 分别代表 A, B, C, p
    // c1, c2, c3 分别代表 x1, x2, x3
    LL k1 = ((c2 - c1 + mod2) * qpow(mod1, mod2 - 2, mod2)) % mod2;
    LL x4 = c1 + k1 * mod1;
    LL k4 = ((((c3 - x4) % mod3 + mod3) % mod3) * qpow((mod1 * mod2) % mod3, mod3 - 2, mod3)) % mod3;
    return (x4 + k4 * ((mod1 * mod2) % mod)) % mod;
}
int main() {
	int n, m, p = 0;
	scanf("%d%d%lld", &n, &m, &mod);
	for (int i = 0; i <= n; i++) scanf("%lld", A + i), A1[i] = A2[i] = A3[i] = A[i];
	for (int i = 0; i <= m; i++) scanf("%lld", B + i), B1[i] = B2[i] = B3[i] = B[i];
	while ((1 << p) <= n + m) p++;
    for (int k = 1; k < 22; k++) {
        for (int i = 1, End = (1 << k); i < End; i++) {
            rev[k][i] = (rev[k][i >> 1] >> 1) | ((i & 1) << (k - 1));
        }
    }
	NTT(A1, rev[p], 1 << p, 1, mod1); NTT(A2, rev[p], 1 << p, 1, mod2); NTT(A3, rev[p], 1 << p, 1, mod3);
	NTT(B1, rev[p], 1 << p, 1, mod1); NTT(B2, rev[p], 1 << p, 1, mod2); NTT(B3, rev[p], 1 << p, 1, mod3);
    for (int i = 0, End = (1 << p); i < End; i++) {
        C1[i] = A1[i] * B1[i]; C2[i] = A2[i] * B2[i]; C3[i] = A3[i] * B3[i];
    }
	NTT(C1, rev[p], 1 << p, -1, mod1); NTT(C2, rev[p], 1 << p, -1, mod2); NTT(C3, rev[p], 1 << p, -1, mod3);
    for (int i = 0; i <= n + m; i++) {
        printf("%lld%c", F(C1[i], C2[i], C3[i], mod), " \n"[i == n + m]);
    }
	return 0;
}

拆系数 FFT

令 $k = \sqrt{p}$（这里可以固定为 32768） ，多项式 $f(x)=k \cdot A(x)+B(x)$，$g(x)=k \cdot C(x)+D(x)$，然后 $f(x) \cdot g(x) = k^2 \cdot A(x)C(x) + k \cdot (A(x)D(x) + B(x)C(x)) + B(x)D(x)$。需进行 8 次 FFT（4 次将 $A(x), B(x), C(x), D(x)$ 用点值表示法表示，4 次将 $A(x)C(x), A(x)D(c), B(x)C(x), B(x)D(x)$ 从点值表示法转换为系数表示法）。常数很大，于是考虑优化。

为啥要拆系数？是为了将系数的位数变小。原本卷积后产生的最大的数约为 $10^9 \times 10^9 \times 10^5 = 10^{23}$，拆系数后，卷积后产生的最大的数约为 $2^{15} \times 2^{15} \times 10^5 = 10^{14}$。可以使用常规数据类型（ 如 long double ）进行储存。

FFT “二合一”优化

假如要对多项式 $F(x)$ 和 $G(x)$ 进行 DFT，则考虑构造如下多项式 $P(x)$，满足 $P(x) = F(x) + G(x) \cdot \text{i}$ 。

可以发现，$P(\omega_n^k)$ 的第 $x$ 项为

$\begin{aligned} (F_x + G_x \cdot \text{i}) \cdot \omega_n^{xk} & = (F_x + G_x \cdot \text{i}) \cdot (\cos(\frac{2\pi x k}{n}) + \sin(\frac{2\pi x k}{n}) \cdot \text{i}) \\ & = F_x \cdot (\cos(\frac{2\pi x k}{n}) + \sin(\frac{2\pi x k}{n}) \cdot \text{i}) + G_x \cdot (\cos(\frac{2\pi x k}{n}) \cdot \text{i} - \sin(\frac{2\pi x k}{n})) \\ & = (F_x \cdot \cos(\frac{2\pi x k}{n}) - G_x \cdot \sin(\frac{2\pi x k}{n})) + (F_x \cdot \sin(\frac{2\pi x k}{n}) + G_x \cdot \cos(\frac{2\pi x k}{n})) \cdot \text{i} \end{aligned}$

$P(\omega_n^{-k})$ 的第 $x$ 项为

$\begin{aligned} (F_x + G_x \cdot \text{i}) \cdot \omega_n^{-xk} & = (F_x + G_x \cdot \text{i}) \cdot (\cos(\frac{2\pi \cdot (-x k)}{n}) + \sin(\frac{2\pi \cdot (-x k)}{n}) \cdot \text{i}) \\ & = F_x \cdot (\cos(\frac{2\pi x k}{n}) - \sin(\frac{2\pi x k}{n}) \cdot \text{i}) + G_x \cdot (\cos(\frac{2\pi x k}{n}) \cdot \text{i} + \sin(\frac{2\pi x k}{n})) \\ & = (F_x \cdot \cos(\frac{2\pi x k}{n}) + G_x \cdot \sin(\frac{2\pi x k}{n})) + (-F_x \cdot \sin(\frac{2\pi x k}{n}) + G_x \cdot \cos(\frac{2\pi x k}{n})) \cdot \text{i} \end{aligned}$

而 $F(\omega_n^k)$ 的第 $x$ 项为 $F_x \cdot \cos(\frac{2\pi x k}{n}) + F_x \cdot \sin(\frac{2\pi x k}{n}) \cdot \text{i}$ ，$G(\omega_n^k)$ 的第 $x$ 项为 $G_x \cdot \cos(\frac{2\pi x k}{n}) + G_x \cdot \sin(\frac{2\pi x k}{n}) \cdot \text{i}$ 。

所以 $\begin{aligned} F(\omega_n^k) = \frac{1}{2} \cdot (\operatorname{Re}(P(\omega_n^k)) + \operatorname{Re}(P(\omega_n^{-k}))) + \frac{1}{2} \cdot (\operatorname{Im}(P(\omega_n^k)) - \operatorname{Im}(P(\omega_n^{-k}))) \cdot \text{i} \end{aligned}$ 。

以及有

$\begin{aligned} G(\omega_n^k) & = (\frac{1}{2} \cdot (\operatorname{Re}(P(\omega_n^k)) - \operatorname{Re}(P(\omega_n^{-k}))) + \frac{1}{2} \cdot (\operatorname{Im}(P(\omega_n^k)) + \operatorname{Im}(P(\omega_n^{-k}))) \cdot \text{i}) \cdot (- \text{i}) \\ & = \frac{1}{2} \cdot (\operatorname{Im}(P(\omega_n^k)) + \operatorname{Im}(P(\omega_n^{-k}))) - \frac{1}{2} \cdot (\operatorname{Re}(P(\omega_n^k)) - \operatorname{Re}(P(\omega_n^{-k}))) \cdot \text{i} \end{aligned}$

因此，对 $P$ 进行一次 DFT 即可。

这样子，我们就只需要进行一次 FFT 就可以对两个多项式进行 DFT 了。

使用上面的方法，就可以使用 2 次 FFT 求出 $A(x), B(x), C(x), D(x)$ 的点值表达式了。

同样，假如要对多项式 $F(x)$ 和 $G(x)$ 进行 IDFT，则同样可以考虑构造如下多项式 $P(x)$，满足 $P(x) = F(x) + G(x) \cdot \text{i}$ 。

对 $P$ 进行一次 IDFT 即可。这样，多项式 $F(x)$ 就在实部，$G(x)$ 就在虚部。

这样子，我们就只需要进行一次 IFFT 就可以对两个多项式进行 IDFT 了。

使用上面的方法，就可以使用 2 次 FFT 求出 $A(x) \cdot C(x),\ A(x) \cdot D(x),\ B(x) \cdot C(x),\ B(x) \cdot D(x)$ 的点值表达式了。

总共只用了 4 次 FFT。而通常来讲，2 次 FFT 和 3 次 NTT 时间差不多（亲测），因此 4 次 FFT 比 9 次 NTT 快出 3 个 NTT（或者 2 个 FFT），比较快。

不过 FFT 对精度的要求较高，因此不用加优化的地方不用乱加，防止精度出问题。

参考代码：（因为精度问题而被卡成了 $50\ pts$）

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1 << 20 | 1;
struct Complex {
	long double x, y; // x + yi
	Complex() {x = y = 0.0;}
	Complex(long double xx, long double yy) {x = xx; y = yy;}
};
Complex operator + (Complex A, Complex B) {
	return Complex(A.x + B.x, A.y + B.y);
}
Complex operator - (Complex A, Complex B) {
	return Complex(A.x - B.x, A.y - B.y);
}
Complex operator * (Complex A, Complex B) {
	return Complex(A.x * B.x - A.y * B.y, A.x * B.y + A.y * B.x);
}
LL t1[N], t2[N], mod, t3[N];
Complex A[N], B[N], C[N], D[N];
Complex AC[N], AD[N], BC[N], BD[N];
int Rev[N];
void FFT(Complex *F, int n, int Sign) { // Sign = 1: FFT , Sign = -1: IFFT
	for (int i = 0; i < n; i++) {
		if (i < Rev[i]) swap(F[i], F[Rev[i]]);
	}
	for (int len = 2; len <= n; len <<= 1) {
		int mid = len >> 1;
		Complex wl = Complex(cos(M_PI / mid), Sign * sin(M_PI / mid));
		for (int i = 0; i < n; i += len) {
			Complex Ang = Complex(1, 0);
			for (int j = 0; j < mid; j++) {
				Complex tx = F[i + j], ty = F[i + j + mid];
				F[i + j] = tx + Ang * ty; F[i + j + mid] = tx - Ang * ty;
				Ang = Ang * wl;
			}
		}
	}
    if (Sign == -1) {
        for (int i = 0; i < n; i++) F[i] = Complex(F[i].x / n, F[i].y / n);
    }
}
Complex P[N];
void PairFFT(Complex *F, Complex *G, int n) {
    for (int i = 0; i < n; i++) P[i] = Complex(F[i].x, G[i].x);
    FFT(P, n, 1); P[n] = P[0];
    for (int i = 0; i < n; i++) F[i] = Complex(0.5 * (P[i].x + P[n - i].x), 0.5 * (P[i].y - P[n - i].y));
    for (int i = 0; i < n; i++) G[i] = Complex(0.5 * (P[i].y + P[n - i].y), -0.5 * (P[i].x - P[n - i].x));
}
void PairIFFT(Complex *F, Complex *G, int n) {
    for (int i = 0; i < n; i++) P[i] = F[i] + G[i] * Complex(0, 1);
    FFT(P, n, -1);
    for (int i = 0; i < n; i++) F[i] = Complex(P[i].x, 0), G[i] = Complex(P[i].y, 0);
}
int main() {
	int n, m, p = 0;
	scanf("%d%d%lld", &n, &m, &mod);
	for (int i = 1; i <= n + 1; i++) scanf("%lld", t1 + i);
	for (int i = 1; i <= m + 1; i++) scanf("%lld", t2 + i);
	while ((1 << p) <= n + m) p++;
	for (int i = 1, End = (1 << p); i < End; i++) {
		Rev[i] = (Rev[i >> 1] >> 1) | ((i & 1) << (p - 1));
	}
    for (int i = 1; i <= n + 1; i++) A[i - 1] = Complex(t1[i] >> 15, 0), B[i - 1] = Complex(t1[i] & 32767, 0);
    for (int i = 1; i <= m + 1; i++) C[i - 1] = Complex(t2[i] >> 15, 0), D[i - 1] = Complex(t2[i] & 32767, 0);
    PairFFT(A, B, 1 << p); PairFFT(C, D, 1 << p);
	for (int i = 0, End = (1 << p); i < End; i++) {
        AC[i] = A[i] * C[i]; AD[i] = A[i] * D[i];
        BC[i] = B[i] * C[i]; BD[i] = B[i] * D[i];
	}
    PairIFFT(AC, AD, 1 << p); PairIFFT(BC, BD, 1 << p);
	for (int i = 0; i <= n + m; i++) {
		printf("%lld%c", (((((1ll << 30) * ((LL)(AC[i].x + 0.5) % mod)) % mod) + ((1ll << 15) * (((LL)(AD[i].x + 0.5) + (LL)(BC[i].x + 0.5)) % mod) % mod) + (LL)(BD[i].x + 0.5)) % mod + mod) % mod, " \n"[i == n + m]);
	}
	return 0;
}

【模板】分治 FFT

题目链接

题目描述

给定序列 $g_{1, \cdots n - 1}$ ，求序列 $f_{0, \cdots n - 1}$ ，满足 $f_i = \sum_{j = 1}^i f_{i - j} g_j$，其中 $1 \le i \le n - 1$，$f_0 = 1$ 。答案对 $998244353$ 取模。

数据范围 $1 \le n \le 10^5$，$0 \le g_i < 998244353$ 。时限 $5s$ 。

第一个整数表示 $n$ ，后面 $n - 1$ 个整数，表示 $g_1, g_2, \cdots g_{n - 1}$。

样例1输入：4\n3 1 2 ，输出：1 3 10 35

题目解答

考虑分治。发现题目的要求类似于卷积，于是考虑使用 FFT/NTT 。

但是后面的数字基于前面的数字，无法快速计算，时间复杂度退化至 $O(n^2)$。

于是我们考虑将类似的转移同时进行，来节省复杂度。

考虑利用分治。

（以下默认 $n = 2^m, m \in \mathbb{N}$ 。对于 $n$ 非 $2$ 的幂次的情况，在序列后面补 $0$ 即可）

在算一段区间 $[l, r)$ 时，将其拆成两个长度一样的区间 $[l, mid)$ 和 $[mid, r)$ ，其中 $mid = \frac{l + r}2$ 。先求出 $[l, mid)$ 之间的答案，之后去计算 $[l, mid)$ 内的 $f$ 对 $[mid, r)$ 内的 $f$ 的贡献。假设其对 $f_x$ 的贡献为 $w_x$ ，其中 $x \in [mid, r)$ ，则有：$w_x = \sum_{i \in [l, mid)} g_{x - i} \cdot f_i$ 。这部分可以利用卷积来快速计算。计算完以后，答案直接加到答案数组就可以了。

我们可以用样例来方便理解：给定 $g = $ [0, 3, 1, 2] ，求 $f$ 。（这里 $g_0 = 0$）

一开始，$f_0 = 1$ 。$f = $ [1, 0, 0, 0]

将其分为两段：[1, 0|0, 0]

先算左边的：[1, 0] ，将其分开：[1|0] 。因为左半边的长度为 $1$，于是不往下递归。计算左区间对右区间的贡献。就是将 [1, 0] 和 $g$ 的前 $2$ 项 [0, 3] 做卷积，得到 [*, 3] ，用 * 表示这个数我们不在乎。将卷积的后半段 [3] 加到这个区间的右半边即可。

操作后得到：[1, 3] ，因为右半边长度为 $1$，于是不往下递归，这区间就做完了，回到上一步。

现在 $f = $ [1, 3|0, 0] ，和上面一样，计算左区间对右区间的贡献。就是将 [1, 3, 0, 0] 和 $g$ 的前 $2$ 项 [0, 3, 1, 2] 做卷积，得到 [*, *, 10, 5] ，将卷积的后半段 [10, 5] 加到这个区间的右半边即可。

操作后得到：[1, 3, 10, 5] ，现在开始计算这个区间的右半段 [10, 5] 。

将其分开：[10|5] 。和之前一样，因为左半边的长度为 $1$，于是不往下递归。计算左区间对右区间的贡献。因为是左区间对右区间的贡献，所以在算卷积时只能要前半段的 $f$ ，后半段均设为 $0$。也就是说，贡献就是将 [10, 0] 和 $g$ 的前 $2$ 项 [0, 3] 做卷积，得到 [*, 30] ，将卷积的后半段 [30] 加到这个区间的右半边即可。

操作后得到：[10|35] 因为右半边长度为 $1$，于是不往下递归，这区间就做完了，回到上一步。

因此得到最终答案：[1, 3, 10, 35] 。

这里贴一下代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N = 1 << 22 | 5;
const LL mod = 998244353ll;
LL A[N], B[N]; // g, f
LL fA[N], fB[N], fAns[N];
int rev[23][N];
inline LL qpow(LL n, LL p, LL mod = 998244353ll) {
    LL ret = 1ll;
    while (p) {
        if (p & 1) (ret *= n) %= mod;
        (n *= n) %= mod; p >>= 1;
    }
    return ret;
}
void NTT(LL *F, int *Rev, int n, int Sign) { // Sign = 1: NTT , Sign = -1: INTT
	for (int i = 0; i < n; i++) {
		if (i < Rev[i]) swap(F[i], F[Rev[i]]);
	}
	for (int len = 2; len <= n; len <<= 1) {
		int mid = len >> 1;
        LL wl = qpow(3, (mod - 1) / len);
		for (int i = 0; i < n; i += len) {
			LL Ang = 1;
			for (int j = 0; j < mid; j++) {
				LL tx = F[i + j], ty = (F[i + j + mid] * Ang) % mod;
				F[i + j] = (tx + ty) % mod; F[i + j + mid] = (tx - ty + mod) % mod;
				Ang = (Ang * wl) % mod;
			}
		}
	}
    if (Sign == -1) {
        for (int i = 1; i < (n >> 1); i++) swap(F[i], F[n - i]);
        LL Invn = qpow(n, mod - 2);
        for (int i = 0; i < n; i++) {
            (F[i] *= Invn) %= mod;
        }
    }
}
void solve(int l, int r, int Log2) { // [l, r)
    if (Log2 == 0) return;
    int mid = l + ((r - l) >> 1);
    solve(l, mid, Log2 - 1);
    for (int i = l; i < r; i++) fA[i - l] = A[i - l];
    for (int i = l; i < mid; i++) fB[i - l] = B[i];
    for (int i = mid; i < r; i++) fB[i - l] = 0;
    for (int i = (1 << (Log2)); i < (1 << (Log2 + 1)); i++) fA[i] = fB[i] = 0;
	NTT(fA, rev[Log2 + 1], 1 << (Log2 + 1), 1); NTT(fB, rev[Log2 + 1], 1 << (Log2 + 1), 1);
    for (int i = 0, End = (1 << (Log2 + 1)); i < End; i++) {
        fAns[i] = (fA[i] * fB[i]) % mod;
    }
    NTT(fAns, rev[Log2 + 1], 1 << (Log2 + 1), -1);
    for (int i = mid; i < r; i++) (B[i] += fAns[i - l]) %= mod;
    solve(mid, r, Log2 - 1);
}
int main() {
	int n, m, p = 0;
	scanf("%d", &n);
	for (int i = 1; i < n; i++) scanf("%lld", A + i); B[0] = 1;
	while ((1 << p) <= n + n) p++;
    for (int k = 1; k < 22; k++) {
        for (int i = 1, End = (1 << k); i < End; i++) {
            rev[k][i] = (rev[k][i >> 1] >> 1) | ((i & 1) << (k - 1));
        }
    }
    solve(0, 1 << p, p);
    for (int i = 0; i < n; i++) {
        printf("%lld%c", B[i], " \n"[i == n - 1]);
    }
	return 0;
}